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INTRODUCTION 

The classical solutions of the two-body problem separate naturally into the 
three cases of elliptic, parabolic, and hyperbolic motion, the mathematics being 
considerably different for each case. A unified formulation is possible, valid 
for all three cases, if certain transcendental functions, which we call the C and S 
functions, are introduced. 

The unified formulation is fully developed by Battin [l] and will not concern 
us here. The purpose of this paper is the presentation of approximations for the 
C and S functions and their derivative functions which reduce significantly the 
computation times required for their evaluation when compared to that required 
by Taylor series expansions. 


THE C AND S FUNCTIONS 

The C and S functions are defined by 

S(x) = (x 1/2 - sin x 1/2 )/x 3/2 , x > 0 (1) 

= [sinh(-x) 1/2 - (-x) 1/2 ]/(-x) 3/2 , x < 0 (2 ) 

C(x) = (l - cos x 1/2 )/x, x > 0 

(3) 

= [l - cosh (-x) 1/2 ]/x, x < 0 (4) 

Since these functions are indeterminate for x = 0 and present accuracy 
problems when evaluated in the neighborhood of x = 0, it is natural to 
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replace the above forms by the following series, convergent for all values 
of x : 


S(x) 


L 

i = n 


(-X ) 1 

( 2i + 3)! 


( 5 ) 


r , , _ V ' 1 (-*) 1 

C ( x ) - (2i + 2)! 

i = 0 


( 6 ) 


For large values of x, the convergence of these series will be slow. It is 
then convenient to use the following reduction formulas, easily derived from the 
definitions (1) through (4): 


in 

X 

i 

n 

X s 

< 

(?) 

2C(4x) = [A(x>] , 

(8) 

4S(4x) = S(x) + A(x) C(x) . 

( 9 ) 


THE C T AND S' FUNCTIONS 

The derivatives S'(x) and C'(x) are needed for certain problems of orbit 
determination, guidance, and optimization. From (1) through (4) we obtain 

S'(x) = [C( x) - 3S(x)]/2x , 

C'(x) = [A( x) - 2C(x)]/2x . 


These forms suffer accuracy problems in the neighborhood of x - 0, 
again forcing us to series representations. Differentiating (5) and (6), 
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we have 


I i(-x ) 1 

S O) = ^ ( ’ 21 +3)! > ( 10 > 

i = 1 


coo = • 01) 

i = 1 


convergent for all values of x . 

For large values of x, the following reduction formulas (obtained by dif- 
ferentiating (7), (8), and (9)) are useful. 


B(x) = 

S(x) + xS'(x) , 

(12) 

C'(4x) 

= -A(x)B(x) 

(13) 

S'(x) 

+ A(x)C'(x) - B(x)C(x) 

(14) 


THE FIKE-KNUTH ALGORITHM 

Our first step in obtaining economical computing forms for (5), (6), (10), and 
(11) was the construction of sixth degree polynomial approximations on various 
intervals in the sense of Chebyshev. In other words, these polynomials minimize 
the magnitude of the maximum error on the interval. The program to accomplish 
this was written by the third author, based on ideas of Stoer [2]. The coefficients 
of these polynomials are given in Numerical Results. 

Assume the approximating polynomial has the form 

P(x) ~ a 0 + ajX + a 2 x 2 + a 3 x 3 + a 4 x 4 + a s x 5 + a 6 x 6 . (15) 

The evaluation of (15) by the usual method of nested multiplication 
requires 6 multiplications and 6 additions. However, using recently 
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developed polynomial evaluation methods [3, 4 ], (15) can be evaluated with 4 
multiplications and 7 additions. The form and parameters for the algorithm, 
as it applies to our functions, are given in Numerical Results. 

In the following description of the algorithm, a 6 is assumed to be positive. 

If a 6 is negative, a minor change is necessary. 

6 . 

Fike's modification of Knuth's method begins with a conversion: letp. - y a 6 , 


t c k - 

a k //i k for k - 0, 1, . . , 5. 

Then compute 


p = 

7 (c, - 1 ) 

D" 

- c 2 - 

PC' 

B' = 

c 4 - p (p + 1) 

E' 

= 2D' 

- B' + 1 

C' = 

c 3 - PB' 

E" 

= 2D" 

- B'D' 

D' = 

P - B * 

E'" 

= c i “ 

B'D" 


Find a real root q of the cubic equation* 

2q 3 + E* q 2 + E" q + E‘" = 0 (16) 


and compute 


A = -y B ' " *1 

C = p - 2A 

B = q - 2AC - A 2 

D = C' - q (l + D') - q 2 - D" - A 2 ( 1 + C) - BC 

E = q 2 + qD' + D" - (A 2 +b)c 

F = c 0 - (q 2 +qD' + D") [c' - q (l + D') - q 2 - D"] . 


*See Appendix I. 
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Then our polynomial can be evaluated as follows: 


qi = /xx 

q 2 = (q, + A ) 2 

^3 = (q 2 + B )(<il + C ) 

P(x) - (q 2 + q 3 + D) (q 3 + E) + F 


In case a 6 < 0, let T(x) - -P(x) and perform all the steps above, except the 
last, for T(x) . The last step should be 

P(x) = “ [T( x )] = (q 2 + q 3 + d) ( _ q 3 - E) ~F . 

If the machine being used has a "load negative" feature which is equivalent 
in execution time to "load positive", and if subtraction is likewise equivalent to 
addition, then this modification is equivalent to the original. 

As Fike points out, his method is a slight variation of that of Knuth [4] , and 
since Knuth's method was inspired by Motzkin [5], the three types bear a strong 
family resemblance. Each begins with a polynomial in form (15) with a 6 = 1 , 
and solves for the parameters in the final evaluation scheme by expanding the 
scheme into a sixth degree polynomial and equating its coefficients with those 
of form (15). To admit treatment of the general polynomial of degree six, how- 
ever, some transformation must be made so that a 6 = 1 . The most straight- 
forward way is 


Q(x) - P(x)/a 6 

and then applying any of the three methods to Q(x) , adding an extra step at the 
last in multiplying the result by a fi . Fike specifies a different sort of trans- 
formation; his may be thought of as converting form (15) into 

^5 ^4 ^3 ^ 2 

R(x) - x 6 + — - x 5 + 7 x 4 + r x 3 + ~ x 2 + “77 x + a n . 

fJL b fl (1 fj , ^ 
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Again, any of the three methods apply to R(x) and values of P(x) are obtained 
by using p-x in the scheme for R(x) , since R(m x ) ~ P(x). 

This transformation, though a bit more complicated, is admirably suited to 
our particular problem. The type of polynomial with which we are dealing has 
the not uncommon characteristic that 

l a 6 l < I a s I < • • • < |a 0 | 

and, in addition, | a 6 | is very small. For example, suppose the coefficients of 
form (15) are 

a 6 = + . 11 x 10' 10 

a s = - . 21 x 10' 8 

a 4 = + . 28 x 10" 6 

a 3 = - . 25 x 10' 4 

a 2 = +. 14 x 10' 2 

a j = - . 42 x 10' 1 

a 0 = + .50 


If we use the division transformation, the coefficients b 4 of Q(x) are 

b 6 = + 1.0 

b 5 = - . 18 x 10 3 

b 4 = + . 24 x 10 s 

b 3 = - . 22 x 10 7 
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b 2 = + .12 x 10 9 

bj = - .36 x 10 10 

b Q = + .44 x 10 11 


Here the errors in the numbers a i have become greatly magnified; worse yet, 
the arithmetic of parameter production using the large numbers b i is likely to 
suffer the effects of large error propagation. In contrast, Fike's transformation 
gives us 

c 6 ~ 1 


C 5 = 

- . 27 x 

10 

C 4 = 

+ . 54 x 

10 

C 3 = 

- .73 x 

10 

C 2 = 

+ . 62 x 

10 

C 1 = 

- . 28 x 

10 

C 0 

.50 



These numbers of manageable size lend themselves very well to whichever 
scheme we choose. For comparison, the two transformations above were 
evaluated by the Knuth algorithm for 40 points over the interval [-1, + l] , 
and the differences between these values and the true values of the polynomial 
were obtained. For the division transformation, the absolute value of the 
maximum error was .92 x 10“ 12 ; for the Fike transformation, this was 
. 16 x 10" 14 , a reduction by a factor of more than 500. Several other test 
cases were run, with results which apparently verify the conclusion that the 
Fike transformation used on this type of polynomial has a very definite advan- 
tage. There are, of course, other transformations which produce polynomials 
in which a 6 = 1 . In general, one should use the transformation which keeps 
the coefficients of the transformed polynomial as small as possible. 
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NUMERICAL RESULTS 


The four approximation polynomials were generated for each of the intervals 
[-1, +l], [-2, + 2 ], [-4, + 4 ], [-16, +16], converted to form (15) and param- 
eters for the Fike evaluation scheme were obtained. In each case, the values 
given by the final scheme were tested against "true" values of the original 
function for all multiples of .002 in the interval concerned. The "true" values 
came from expanding the power series of the function (1) for enough terms to 
guarantee that the relative error from truncation would be less than 10“ 15 . 

The following tables exhibit, for each of the sixteen functions considered, the 
coefficients for form (15) , the parameters A, B, C, D, E, and F for the Fike 
scheme, and the maximum absolute errors for both methods. For comparison, 
a degree 4 approximation polynomial was evaluated by both methods for the 
functions C(x) and S(x) on the interval [-1, + 1] , and the results are presented 
here also. 
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C(X) = a 0 + a^ x + a 2 x 2 + a 3 x 3 + a 4 x* + a 5 x 5 + a 6 x 6 


on E"^ 1/ 0 


h = 1 

a 0 = +0.4999999999999998 X 10° 
a, = -0.4166666666667176 X 10" 1 
a 2 = +o. 1388888888888999 x 10" 2 
a 3 = -0.2480158725995993 X 10" 4 
a 4 = +0.2755731917059028 X 10" 6 
a 5 = -0.2087759200397967 X 10" 8 
a 6 = +0.1 14713410831 1665 X 10" 10 

M a = 0.76327833 X 10" 15 


h = 4 

a 0 = +0.4999999999998401 X 10° 
a , = -0.4166666668808485 X 10" 1 
a 2 = +0. 1388888889138842 x 10" 2 
a 3 = -0.2480157659289839 X 10" 4 
a 4 = +0.2755731272513377 x 10" 6 
a 5 = -0.2089014196095935 X 10" 8 
a 6 = +0. 1147636934013430 X 10' 10 

M a = 0.12239224 X 10" 10 


h = 2 

a 0 =+0.4999999999999993 X 10° 
a, =-0.4166666666700118 X 10" 1 
a 2 = + 0. 1388888888892785 X 10" 2 
a 3 = -0.2480158663241807 X 10" 4 
a 4 =+0.2755731881710992 X 10' 6 
a 5 = -0.2088010277268315 X 10" 8 
a 6 =+0.1147215380312168 X 10" 10 

M a = 0.95645714 X 10“ 13 


h = 16 

a 0 = +0.4999999894793170 X 10° 
a, = -0.4166675473500692 X 10" 1 
a 2 = +0. 1388889916034137 X 10 -2 
a 3 = - 0.24798836331 19184 x 10" 4 
a 4 = +0.2755565077419917 X 10" 6 
a 5 = -0.2109148028487573 X 10' 8 
a 6 = +0.1156091702389399 x 10" 10 

M a = 0.20159773 X i<r 6 
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on [" h, hj 


9 ] = 



x 


92 = (^1 + A) 2 
^3 = (^2 + B) (9i + C) 
C(X) = (92 + 93 + D) (+ 93 + E) + F 


h= 1 

A = +0.4513408582627891 X 10° 

B = +0.3744865190483202 X 10 1 
C = -0.2769272800754423 x 10 1 
D = + 0.9433565393074166 x 10 1 
E = + 0. 1055413968372178 x 10 2 
F = +0.6288190624578802 x 10" 2 

M = 0 . 5828670879282069 X 1 0" 1 4 


h = 2 

A = + 0.4512008438957284 X 10° 

B = +0.3743936586512442 x 10 1 
C = -0.2769076432349482 x 10 1 
D = + 0.9429681327692907 x 10 1 
E = +0.1055053194443676 x 10 2 
F = +0.6284207351324604 x 10" 2 

M = 0.9935108291614367 x io -13 


h = 4 

A = + 0.4507019590625572 X 10° 

B = +0.3740448131455307 X io 1 
C = -0.2768317205414987 x IO 1 
D = +0.9414899439055362 x 10 1 
E = +0.105370131 1149719 X 10 2 
F = +0.6272339359526764 x io - 2 

M = 0.1 224323420423443 x l o' 1 0 


h = 16 

A = + 0.4405766736959988 x 10? 

B = + 0.3669989101432331 X 10 
C = - 0.2752824996097087 x 10 1 
D = + 0.91 17484138879304 X io 1 
E = +0. 1026464040151929 X 10 2 
F = +0.6161613003116790 x io' 2 

M = 0.2015977389469012 x io -6 
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S( X) = a 0 + aj x + a 2 x 2 + a 3 x 3 + a 4 x 4 + a$x 5 + a 6 x 6 


on j-* h, hj 


h = 1 

a 0 = +0. 1666666666666665 X 10° 
a, = - 0. 8333333333333568 x 10" 2 
a 2 = +0. 1984126984129264 x 10" 3 
a 3 =-0.2755731919939401 x 10"^ 
a 4 = + 0. 2505210785999854 x 10"* 
a 5 = - 0. 1605953765319026 X 10" 9 
a 6 = + 0. 7650283228592385 X 10" 1: 

M a = 0.55511151 X 10" 16 


h = 4 

a 0 = + 0. 1666666666666581 X 10° 
a, =-0.8333333334593103 X 10" 2 
a 2 = +0. 1984126984258407 X 10' 3 
a 3 = - 0.2755731292513204 X 10' 5 
a 4 = +0.2505210496761327 X io -7 
a 5 = -0. 1606691704048320 X 10' 9 
a 6 = + 0.7650122280184766 X 10' 12 

M a = 0.72000739 x 10" 12 


h = 2 

a 0 = + 0. 1666666666666663 X 10° 
a, = -0.8333333333352921 x 10" 2 
a 2 = + 0. 1984126984129057 X IO' 3 
a 3 = -0.2755731883077317 x 10' 5 
a 4 = + 0.2505210816170333 X 10“ 7 
a 5 = - 0.1606101 133600677 X 10 " 9 
a 6 = + 0 . 7647926042737674 X l O' 12 

M a = 0.56621374 x 10" 14 


h = 16 

a 0 = + 0.1666666661133027 X 10° 
a , = - 0 . 8333338509758059 x l O' 2 
a 2 = + 0. 1984127524406629 X IO' 3 
a 3 = -0.2755570214946503 x 10" 5 
a 4 = + 0.2505123071826789 x 10" 7 
a 5 = -0.1618528418504030 x 10" 9 
a 6 = + 0.7694603615375217 X 10' 12 

M a = 0.11845921 X 10" 7 
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on 


q 2 = (q, +a) 2 
^3 = (^2 + b) (q, +c) 

S(X) = (q 2 +q 3 + D ) ( + q 3 +E) + F 



h = 1 

A = + 0. 1030541 1 10544949 x 10° 

B =+0.1357446199107850 X 10 1 
C = -0. 1709888012144409 x 10 1 
D = +0. 1803960616654583 x 10 1 
E =+0.2062171852872241 xio 1 
F = +0.2130010466488949 x 10" 1 

M = 0 . 2359223927328455 x 1 0“ 1 4 


h = 2 

A = + 0. 1028274494783523 x 10° 

B = +0. 1356879505417743 x 10 1 
C = - 0. 1709784631095070 x 10 1 
D = + 0. 1802832347589989 x 10 1 
E = + 0.2060949727430426 X 10 1 
F = + 0 . 2128753997322974 x 10" 1 

M = 0.8076872504148012 x 10" 14 


h = 4 

A = + 0. 1026453527945748 x 10° 

B = + 0.1356011711587295 x 10 1 
C = -0.1709549340846876 X 10 1 
D =+0.1800557277079013 x 10 1 
E = + 0.2059246874696585 X 10 1 
F = + 0.2125209137884825 x 10" 1 

M = 0.7223804887601657 XIO' 1 


h= 16 

A = + 0. 9898746283297480 x 10“ 1 
B = + 0. 1338586121335949 x 10 1 
C = -0. 1704756186376375 X 10 1 
D = + 0. 1754906650979839 x 10 1 
E = + 0. 2025050653157090 x 10 1 
F = + 0. 2056593593968585 x 10" 1 

M = 0.1 184592103575797 x 10“ 7 
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C'( X) = a 0 + aj x + a 2 X 2 + a 3 x 3 + a 4 x 4 + a 5 x 5 + a 6 x 6 on [- h, hj 


h = 1 

a 0 = - 0.4166666666666430 X 10" 1 
a, = + 0.2777777778467318 X 10“ 2 
a 2 = - 0. 1488095238678453 X 10' 3 
a 3 = + 0.6613751098214413 X 10' 5 
a 4 =-0.2505208412532178 x 10' 6 
a 5 = + 0.8269965205281566 x 10" 8 
a 6 = - 0.2412214891589441 x 10' 9 

M a = 0.21684043 X10" 16 


h = 4 

a 0 = - 0.4166666651 125364 X 10" 1 
a, = +0.2777780643726382 X 10' 2 
a 2 = -0. 1488097663598733 x 10" 3 
a 3 = +0.6612326049924104 x 10' 5 
a 4 = -0.2504581368842600 X 10' 6 
a 5 = +0.8437138609417991 x 10' 8 
a 6 = -0.2463133603265637 x 10" 9 

M a = 0.31997321 X 10' 12 


h = 2 

a 0 =-0.4166666666606742 X 10" 1 
a, =+0.2777777822034584 X 10' 2 
a 2 = - 0. 1488095275535492 X 10" 3 
a 3 = +0.6613668137463213 X 10" 5 
a =-0.2505171918626241 X 10" 6 
a ^ = + 0 . 8303 1 36595603477 X 10' 8 
aj = -0.2422316681583509 X 10’ 9 

M a = 0.25058081 X 10' 14 


h = 16 

a 0 =-0.4166666641763858 X 10' 1 
a, = + 0.2777780078584500 X 10' 2 
a 2 = -0.7440478621862503 X JO' 4 
a 3 =+0.1102220894089232 X 10' 5 
a 4 = - 0. 1043798352532477 X 10' 7 
a 5 = + 0.6938557071056230 X l O' 10 
a 6 =- 0.3366982823031963 X 10" 12 

M a = 0.52654049 XI O' 8 
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h = 1 


^1 - v^a 6 x 
q 2 = (9i + A) 2 
q 3 =(q 2 + B) (q, +C) 
C'(X) = (q 2 +q 3 + D)(-q 3 -E)-F 



A = +0.6421500675794880 x 10 - 1 
B = +0.9631412054424315 X 10° 

C = - 0. 1485378935681960 X 10 1 
D = +0. 1206516480446427 x 10 1 
E = +0. 1262370275488431 X 10 1 
F = +0.2235785774036898 x 10 " 2 

M =0. 1786765180256109 X 10' 15 


h = 2 

A = + 0.6415596248585610 x lo -1 
B = + 0. 9629599640982438 x 10 0 
C = - 0. 1485350658374242 x 10 1 
D = +0. 1206153148806730 X 10 1 
E = +0. 1262088821410514 X 10 1 
F = +0.2230745797490601 X 10 " 2 

M = 0.2680147770384164 x 10" U 


h =4 

A = +0.6405925631551870 x 10" 1 
B = +0.9624584033933687 x 10 ° 
C = -0.1485269967153406 X 10 1 
D = +0.1205018912730453 x 10 1 
E =+ 0 . 1261394717119364 x 10 1 
F = +0.2210887038704210 x 10" 2 

M = 0.3200469403386028 x 10' 1 


h = 16 

A = + 0.6208183431484613 x 10- 1 
B = + 0.9523193847619871 x 10° 

C = -0.1483582934430130 x 10 1 
D = + 0 . 1 182131 124018069 x 10 1 
E = + 0.1247277466997771 x 10 1 
F = + 0. 1829570481293727 x 10 “ 2 

M = 0.5265405070287163 x 10" 8 
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S'(X) = a 0 + a, X + a 2 x 2 + a 3 x 3 + a 4 x 4 + a 5 x 5 + a 6 x 6 


on £- h, hj 


h = 1 

a 0 = -0.8333333333333210 X 10 * 2 

a, =+0.3968253968616768 x 10 * 3 

a 2 = -0.1653439153716376 X 10~ 4 

a 3 = +0.60125031 10063301 X 10* 6 

a. = - 0 . 1927084107177350 x 10" 7 

a 5 = + 0.551 1761621292874 x 10' 9 

a . = -0.1418571972378231 x 10* 10 
6 

M a = 0.26020852 X 10* 17 


h = 4 

a 0 = -0.8333333325953303 X 10 * 2 
a, = +0.3968255472561213 x 10 * 3 
a 2 = -0. 1653440305444520 X 10 * 4 
a 3 = +0.601 1754909568150 X 10 * 6 
a 4 = -0. 1926786201990815 x 10 * 7 
a 5 = +0.5599576811597955 X 10 * 9 
a 6 = -0.1442776137237285 X lO' 10 
M a = 0.16841563 x 10" 13 


h =2 

a 0 = -0.8333333333304809 X 10* 2 
a 1 = +0.3968253991531185 x 10* 3 
a 2 = -0. 1653439171256403 X 10* 4 
a 3 = + 0.6012459474285134 x 10* 6 
a 4 = -0.1927066737597982 X 10* 7 
a 5 =+0.5529210296030711 x 10 * 9 
a 6 = -0. 1423381311296310 X 10* 10 

M a = 0.13444107 X 10* 15 


h = 16 

a 0 = -0.8333333321480760 x 10* 2 
a. = + 0.3968255178485000 X 10" 3 
a 2 = -0.8267196924462379 X 10 -5 
a 3 = + 0. 1002046526767437 X 10' 6 
a 4 = -0.8029333915166488 X 10" 9 
a 5 = +0.4617817733427159 X 10* 11 
a 6 = -0.1978183008506138 x 10" 13 

M a = 0.27690121 x 10" 9 
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9 i = VaJ x 

q 2 = (^i A ) 2 

q 3 = (q 2 + B) ( qi + C) 

S'(X) = (q 2 + ^3 + D ) (“^3 “ E) _ F on J- h, hj 


h = I 

A =-0.4272502910304863 X|0 _1 
B =+0.4460409572215307 X |0° 

C =-0.1020276052465536 X I0 1 
D =+0.4064530281034859 X 10° 

E =+0.3450015283629582 X 10° 

F =+0.2885054406289098 x I0' 2 

M= 0.5290906601729259 X|0" 16 


h = 4 

A =-0.4274056995233864 X |0" 1 
B =+0.4458621272744412 x 10° 

C =-0.1020307887595762 x I0 1 
D =+0.4061739766937458 x 10° 

E =+0.3449170623874602 x 10° 

F =+0.2876661053629565 x |0' 2 

M= 0. 1689186984732414 X|0“ 13 


h = 2 

A =-0.4272336986330691 x|o -1 
B =+0.4460090638507059X|0° 

C =-0.1020282840851638 X 10 1 
D =+0.4064007086884238 X |0° 

E =+0.3449898120978949 X |0° 

F =+0.2883373149053849 X |0' 2 

M= 0.1 821 459649775645 X |0 -15 


h = 16 

A = -0 . 431 0245077956 152 x 1 0“ 1 
B =+0.4428994350437274 x 10° 

C =-0.1020789569928175 X I0 1 
D =+0.4016042135773923 X 10° 

E =+0.3434241762739658 x 10° 

F =+0.2744481265299470 x 10 -2 

M= 0.2769012918263367 x I0‘ 9 
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P(X) = a 0 + a]X+a 2 x 2 +a 3 x 3 +a 4 x 4 


C(X)#4 [-l, + l] 

a 0 = + 0.5000000000007167 X 100 
a, = -0.4166666601424471 X 10' 1 
a 2 = + 0. 1388888879568642 X 10" 2 
a 3 = -0.2480419696201834 X 10" 4 
a 4 = + 0.2755932664579228 X 10' 6 

M 0 = 0.13048673 X 10“ 9 


S(X)*4 [-1, + l] 

a 0 = + 0.1666666666667142 X 10° 
a, = -0.8333333283147393 X 10" 2 
a 2 = + 0. 1984126977913694 X 10" 3 
a 3 = -0.2755932664326337 X 10" 5 
a 4 = + 0.2505344666229896 X 10" 7 

M a = 0. 10037290 X 10“ 10 


q,= A x 

q 2 = K + B ) 2 

p (x) = (q, + q 2 + C) (q 2 + D) + E 


C(X)*4 [-1,+l] 

A = + 0.2291221893900772 x 10” 1 
B =-0.7655416156428518X 10° 
C = + 0 . 2592589041 826887 X 1 0° 
D = + 0.4011 554686880556 X 1 0° 

E = - 0 . 33450083938941 42 X 1 0° 

M = 0.1304868574303339 X 10" 9 


S(X)*4 [-1, +l] 

A = + 0.12581 04947765253 X 1 0" 1 
B = -0.5959855810891701 X 10° 

C = + 0.1104585254341674X10° 

D = + 0.2038526154314646 X10° 

E = - 0 . 9365973340739070 X 1 O' 1 

M = 0.1003741534333355 X IQ" 10 
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APPENDIX I 


The cubic equation (16) was solved by an interval-halving technique which can 
be extended to any continuous function f . In general, this technique gives us the 
"smallest” (in the sense of representability by computer) interval in which a 
value x can lie such that f(x) = 0, and this smallest interval can be found by 
a finite, fixed number of iterations. The most familiar interval-halving process 
consists of two parts: (1) finding the initial bracketing interval, the interval 
[a, b] in which f(x) changes sign, then (2) successively halving this interval, 
choosing each time the subinterval in which f(x) changes sign, until the interval 
is as small as desired. If part (1) is performed properly, then part (2) can be 
performed with a number of iterations determinable a priori, thus eliminating the 
test for interval size at the end of each iteration. 

To begin, choose a number v and a number t, of the same sign as v, such 
that 1 1 1 < | v| . [Positive v's are used for positive roots and negative v's for 
negative roots.] Test the following sequence of intervals for a sign change of 
f(x): 


(1) [v, 

v + t] 



(2) [v + 

t, v + t 

+ 2t] 


(3) [v + 

t + 2t, \ 

+ t + 2t 

+ 4t] 

• 

P 


p+ 1 

(P + 2) 


24, v 

+ E 2it 


i = 0 


i = 0 


until the initial bracketing interval is found. Since 1 1| < | v| , we have 

| v | + 2 q+ 1 | 1 1 - | 1 1 > 2 q+1 1 1 1 => 

+ 1 1 1 (2 q+1 - l) > 2 q+1 1 1 1 => 


v 


v + t 


q 

i =0 


2 q +i t 
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and since v and t have the same sign, we have 



> | 2 q+ 1 t| 


This simply says that the length of the initial bracketing interval is less than, or 
equal to, the magnitude of the small end; in turn, this means that the large end 
of the interval is at most twice the magnitude of the small end. 

Now, consider how the endpoint values would be represented in floating- 
point binary arithmetic (normalized) with an r-bit fraction. If the difference 
between their binary exponents is at most 1 (which is what we are getting at 
above), then it can be seen that the number of distinct points in the initial 
bracketing interval is at most 2 r . Therefore, the number of interval-halving 
iterations needed— that is, the number of times one reduces his choice of points 
in the interval by one-half— is at most r. Moreover, it often turns out that f is 
nearly (or exactly) zero at an end point of one of the half-intervals, so that r 
iterations are not always needed. 

We have treated the special case 1 1| < | v| , but we need not restrict our- 
selves to it. The number of interval-halving iterations needed depends upon the 
size of t, and if one is willing to iterate a bit more he can find the initial 
bracketing interval more quickly by increasing t; the converse of this also 
holds. 
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